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Abstract. 

QJ ' We report a numerical simulation of the rate of crystal nucleation of sodium chloride 

from its melt at moderate supercooling. In this regime nucleation is too slow to be 
studied with "brute-force" Molecular Dynamics simulations. The melting temperature 



Ph. 



C/3 

O ■ of ("Tosi-Fumi") NaCl is ~ 1060K. We studied crystal nucleation at r=800K and 

. 825K. We observe that the critical nucleus formed during the nucleation process has 

the crystal structure of bulk NaCl. Interestingly, the critical nucleus is clearly faceted: 
the nuclei have a cubical shape. We have computed the crystal-nucleation rate using 
two completely different approaches, one based on an estimate of the rate of diffusive 
crossing of the nucleation barrier, the other based on the Forward Flux Sampling and 
^ . Transition Interface Sampling (FFS-TIS) methods. We find that the two methods 

' yield the same result to within an order of magnitude. However, when we compare 

the extrapolated simulation data with the only available experimental results for NaCl 
ly-^ , nucleation, we observe a discrepancy of nearly 5 orders of magnitude. We discuss the 

' possible causes for this discrepancy. 

o 

o\ 
o 

1. Introduction 

•rH 

X 

^ ' Crystallization of salts is a phenomenon of great practical relevance. In fact, it is one 
of the most important industrial separation processes. But it also plays a crucial role in 
geological processes that occur on an altogether different time scale. The crystallization 
process consists of two steps: nucleation and growth. If nucleation is slow compared to 
the time it takes a crystal to grow to a size comparable to the size of the container, large 
single crystals will form (an example is rock salt). When nucleation is fast, the resulting 
solid will form as a fine powder. It is clearly important to be able to predict the rate 
of nucleation of salts and - at a later stage - to understand the factors that influence 
nucleation. In the present paper we aim to demonstrate that, with current simulation 
techniques and currently available force-fields, it is indeed possible to compute the rate 
of nucleation of a real salt crystal (in the present case NaCl from its melt) This opens 
the way to "ab-initio" predictions of nucleation rates of many ionic substances. 
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Solutions or melts can often be cooled well below their freezing temperature. The 
reason is that the formation of small nuclei of the stable crystal phase is an activated 
process that may be extremely slow. Intuitively, it is easy to understand why crystal 
nucleation is an activated process, i.e. why there is a free-energy barrier separating 
the metastable parent phase (the liquid) from the stable crystal phase. The point is 
that, initially, the formation of small crystalline nuclei costs free energy. But once the 
crystal nucleus exceeds a critical size, its free energy decreases as it grows. The rate at 
which crystal nuclei form depends strongly on AGcrit, the free-energy required to form a 
critical nucleus. Classical Nucleation Theory (CNT) is commonly used to estimate the 
height of the nucleation barrier and to predict the rate of crystal nucleation [U El El HI E] . 
According to CNT, the total free energy of a crystallite that forms in a supersaturated 
solution or melt contains two terms: the first is a "bulk" term that expresses the fact 
that the solid is more stable than the supersaturated fluid - this term is negative and 
proportional to the volume of the crystallite. The second is a "surface" term that 
takes into account the free-energy cost of creating a solid-liquid interface. This term is 
positive and proportional to the surface area of the crystallite. According to CNT, the 
total (Gibbs) free-energy cost to form a spherical crystallite with radius R is 

AG = —R''psAfi + A7rR'-f, (1) 

where ps is the number density of the solid, AyU (< 0) the difference in chemical potential 
of the solid and the liquid, and 7 is the solid-liquid interfacial free energy density. The 
function AG goes through a maximum at i? = 27/(p5'| Ayu|) and the height of the 
nucleation barrier is 

AG,„, = _7V(p^|A/.|)l (2) 

The crystal-nucleation rate per unit volume, /, depends strongly on AGcri*: 

I = Kexp{-AGcrit/kBT). (3) 

Here k is a kinetic prefactor, T is the absolute temperature and kB is Boltzmann's 
constant. The CNT expression for the nucleation rate then becomes 



K exp 



3 J 

Under experimental conditions, nucleation is infrequent on the time scale of typical 
molecular processes. Yet, when it happens, it proceeds rapidly. This makes it difficult 
to study the structure and dynamics of crystal nuclei of atoms or small molecules in 
experiments. In the case of NaCl, the experiments are also complicated by the fact that 
crystallization occurs at high temperatures. This may explain why there is a scarcity of 
experimental data on the nucleation of NaCl. To our knowledge, the only data are those 
of Buckle and Ubbelhode P, [71 [8] from the 1960's. In these experiments, crystallization 
in NaCl micro- droplets was observed visually. As the droplet size (((9(3/i)) and time 
window for the measurement (1-30 seconds - after which the droplets sedimented out 
of view) - were fixed, the nucleation rate could be determined at one temperature only 
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(905 K for NaCl). At this temperature, the nucleation rate was such that, on average, 
one nucleus would form during the observation time ((9(10s)) in a droplet with a volume 
of order 10~^^m^. Hence the experimental nucleation rate per unit volume was 0(10^^) 

— s —1 
m s . 

For experimental nucleation rates of this order of magnitude, brute-force MD 
simulations are out of the question. The average time it would take for nuclei to form 
spontaneously in a system consisting of several thousands of particles is of order of 10^° 
seconds. Clearly, this is beyond the scope of MD simulations. The standard solution to 
circumvent this problem is to perform simulations at much larger undercooling than used 
in the experiments. Huang et al. [9] performed MD simulations of melting and freezing of 
a droplet composed of 216 NaCl ions in vacuum: to this end, they performed temperature 
quenches down to 550K (i.e. approximately half the melting temperature) and found 
nucleation rates of the order of O(10^^)m~^s~^, which is 20 orders of magnitude higher 
than the experimental rate at 905/^. 

Another effort to study nucleation at less severe supercooling was made by Koishi et 
al. [To]. These authors performed an MD simulation of 125000 ions system in vacuum, 
at temperatures of 740K (i.e. approximately 0.7 the melting temperature Tm) and 640K 
(i.e. approximately 0.6 the melting temperature T^). Both free and periodic boundary 
conditions were used. The estimated nucleation rate at TAOK was O(10^^)m~^s~^, which 
is virtually the same value as found by Huang et al. [9] at a much larger supercooling. 
This is surprising because nucleation rates tend to depend very strongly on temperature. 
This suggests that, at least at the lowest temperatures, the barrier for crystal nucleation 
is negligible. More in general, crystal nucleation under extreme supercooling need not 
proceed following the same path as under moderate supercooling 

In what follows, we use the technique of refs. [11], [121 US] based on a combination 
of umbrella sampling [13] (to determine the barrier height) and a dynamical simulation 
(to determine the crossing rate). The computing time required for this scheme does 
not scale exponentially with the nucleation barrier, but it does increase with increasing 
nucleus size. 

We also compute the nucleation rate using an algorithm based on the "forward flux 
sampling" (FFS) [15] and the "Transition Path Samplig" (TIS) techniques [16] and 
compare it with the one obtained using the method previously mentioned. 

In the present work we study homogeneous crystal nucleation in the Tosi-Fumi 
NaCl model for NaCl at two different temperatures, viz. Ti=800K and T2=825K, 
corresponding to 25% and 22% undercooling. For this system, we computed the 
nucleation barrier, examined the structure and shape of the critical nucleus and 
computed the nucleation rate. 
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The Tosi-Fumi rigid-ion interaction potential for NaCl is of the following form [171 HH], 

£/„M=A.e'«"«-''«-^-^^ + !^. (5) 
where the parameters have the values given in table [H 
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[kJ/mol] 
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[A^kJ/mol] 


[A^kJlmol] 


\M 


Na-Na 


25.4435 


3.1546 


101.1719 


48.1771 


2.340 


Na-Cl 


20.3548 


3.1546 


674.4793 


837.0770 


2.755 


Cl-Cl 


15.2661 


3.1546 


6985.6786 


14031.5785 


3.170 



Table 1. Potential parameters for NaCl. 



This pair potential is written as the sum of a Born-Mayer repulsion, two attractive 
van der Waals contributions and a Coulomb interaction term. In our simulations, we 
calculated the Coulomb interactions using the Ewald summations method with a real 
space cutoff of 10 A and a real space damping parameter of 0.25 A^^ . We truncated the 
Van Der Waals part of the potential at 9 A, assuming the g{r) = 1 beyond this cutoff. 

The computed number density of ions in the bulk solid at 800K and 825K at 10^ Pa 
was 0.041y4~^, in agreement with experiment [19]. The density of ions in the supercooled 
liquid at the same temperature and pressure was 0.034^4"^. 

We prepared under cubic boundary conditions a supercooled system of (12)'^ NaCl 
ion pairs at ambient pressure by cooling it down below the melting temperature. For the 
present model, Anwar et al. [20] have computed the melting temperature: Tm=(1064± 
14)K, which is very close to the experimental melting temperature (T^^p=1072K). 
Using constant-pressure Monte Carlo simulations, we cooled the system down to the 
temperatures where we studied nucleation: Ti=800K and T2=825K, corresponding to 
25% and 22% undercooling. Note that the experiments on NaCl nucleation where 
performed at a somewhat higher temperature (16% supercooling). The reason why we 
could not perform simulations at these higher temperatures is that the critical nucleus 
would be about twice the size of the nucleus that could be studied without spurious 
finite-size effects [21] for the system sizes that we employed. At temperatures below 
750K, spontaneous nucleation occurred during the simulations. We therefore kept the 
temperature above this lower limit. 

Nucleation is an activated process. In steady-state, the nucleation rate per unit 
volume and time is given by eq. ([3]) [2l [221 E] where exp(— /3AGcrit) is the equilibrium 
probability per nucleus to find a critical nucleus in the metastable parent phase, k is a 
kinetic prefactor. In the case of a diffusive barrier crossing, k can be expressed as: 
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where puq is the ion density of the metastable hquid, f^cru the rate at which particles 
are added to a critical nucleus, An=fiiiq — psoi is the difference in chemical potential 
between liquid and solid and \ ^ — the Zeldovitch factor, that takes into account 
the fact that, during steady-state nucleation, the concentration of critical nuclei is, in 
fact, not the equilibrium concentration. 

It is important to distinguish between fncrity which is the rate at which particles 
are added to a nucleus with the critical size, and the net flux across the nucleation 
barrier. In steady state, this net flux is equal the number of nuclei that go from ncrit to 
ncrit + 1 minus the number that go from ncrit + 1 to ncrit [21 El H]. Hence, the actual 
nucleation rate is a combination of the forward rate and the backward rate /~. 
However, because of detailed balance, knowledge of f ncrit (combined with knowledge of 
the barrier height and shape) is enough to compute the nucleation rate. 

Assuming a diffusive attachment or detachment of single particles from the critical 
nucleus, the forward rate f^^ at the top of the barrier can be related to the spontaneous 
fluctuations in the number of particles in a nucleus at the top of the nucleation barrier: 

J crit 2 t 

where (An^^j^(t)) = {[ncrit{t) — 'n.crit(O)]^) is the mean square change in the number of 
particles belonging to the critical nucleus during a time interval t. To estimate f^^^ a 
series of dynamical trajectories were necessary: after generating a set of uncorrelated 
configurations at the top of the barrier, we carried out NVT MD simulations using the 
DL_POLY package[23j with a timestep of 0.5 fs. We computed the nucleation rate 
using equation Q. 

Moreover, we computed the nucleation rate per unit time and volume using an 
algorithm based on path-sampling techniques [Ul [IS] and compared the results with 
those obtained using equation (j3]). 



3. Results 



We computed the free energy barrier for crystal nucleation at Ti=800K and T2=825K, 
corresponding to /3A/ii=0.54 and /?A/i2=0.48. The values of j3Afj, were estimated 
numerically by thermodynamic integration from the coexistence temperature and free 
energies reported by Anwar et al. [20] . Fig{T] shows the computed nucleation barriers 
as a function of n. 

As expected, AG decreases with supersaturation. Around T<750K the barrier 
gets sufficiently low that spontaneous nucleation can take place on the time scale of a 
simulation. The size of the critical nucleus, N^, was estimated according to a fit of the 
functional form of the CNT, and we found Nc ~ 120 ions at Ti and Nc ~ 150 ions at T2. 
Koishi et al. [10] estimated Nc = 120 — 130 ions at 640K and 740K, which is surprising 
in view of the CNT prediction that the size of the critical nucleus scales as (7/I A/i|)^. If 
we make the usual assumption [5| that 7 ~ Ah and A/i ^ {Ah/Tm){T — Tm), where Ah 
is the entropy of fusion per ion pair, then we would expect that iVc ~ (1 — T /Tm)~^ and 
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Figure 1. Free energy barriers AG as a function of the nucleus size n for 
Ti=800K(/3A//i = 0.54) and r2=825K(/5A/i2 = 0.48). Error-bars on /3AG are of 
the order of 1 ksT. The dashed curves are fits to the functional form given by CNT. 

we would predict that the critical nucleus at 640K should contain a quarter as many 
particles as those at 800K. 

Next, we consider the structure and shape of the critical nucleus. Fig|2] shows a 




Figure 2. Snapshot of the critical nucleus at T2=825K: the bulk NaCl structure is 
already evident. The critical crystal nucleus seems to have rudimentary facets, which 
is in agreement with the interpretation of the experiments in ref . [6] . 

snapshot of the critical nucleus at T2=825K. Note that the crystal presents rudimentary 
low-index facets. In experiments [6], El] the existence of such facets was postulated, as 
they may act as sites for subsequent heterogeneous nucleation. 

As can be seen from FiglJl the critical nucleus already shows the charge-ordered 
rock-salt structure of the bulk phase. This indicates that, at least for NaCl, nucleation 
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does not proceed via an intermediate metastable phase [25]. Figl2] also shows that, in 
the temperature range that we studied, the critical nucleus is non-spherical. In order 
to quantify the degree of non-sphericity of the critical nucleus, we expanded its density 
with respect to the center of mass in rank-four-spherical-harmonics and constructed the 
quadratic invariant 5*4 [26]: we obtained S'4(Ti)=0.115 and S'4(T2)=0.110. For a simple 
cube 5*4=0.172 and for a sphere 5*4=0, therefore the shape of the critical nucleus is closer 
to a cube than to a sphere. In other words, the critical nucleus already exhibits the 
morphology of macroscopic NaCl crystals |27|. 

In order to make sure that there were not finite size effects that resulted in 
interactions between image cluster, we visually checked that the critical nuclei didn't 
show a preferred spacial orientation and that the minimum distance between them was 
bigger than half box. 

Moreover we calculated the Debye-Huckle screening lenght and found that it was 
smaller than 1 A, showing that the critical nuclei did not even electrostatically interact. 
We could therefore conclude that there was no induced nucleation due to the interaction 
between a critical cluster and its own periodic image. 

Using the computed height of the nucleation barrier and the values of A/i as input, 
we can estimate the surface free-energy density 7/s. To this end, we make use of the 
CNT expression for the barrier height (eqn. [2]). However, this expression assumes that 
the critical nucleus is spherical. It is easy to derive the corresponding expression 
for a cubical nucleus. The results for both estimates of are given in Table [2j 
There exist experimental estimates of at 905K [6]. These estimates are based on 



T[K] 


"^sphere 


7cu6e 


800 


98 ±2 


80 ± 1 


825 


99 ± 1 


79 ± 1 


905 


84.l(^) 


67.8('^) 



Table 2. Surface free energy density (in erg cm^'^) assuming spherical and cubical 
shape for critical nuclei. At 905K we report the experimental value 6 . The entry 
in the lower right-hand corner is based on the experimental estimate, but assuming a 
cubical nucleus. 

a somewhat questionable CNT expression for the nucleation rate. Moreover, in ref. [6] 
it is assumed that the critical nucleus is spherical. The experimental estimate of 7/s 
{lexp=84:.l[erg cm^^]) is therefore not based on a direct determination. Nevertheless, 
in the absence of other experimental data, this is the only number we can compare to. 

As the table shows, there is a fair agreement between simulation and experiment. 
The experimental estimate for 7 is based on the assumption that the critical nucleus is 
spherical. If it is cubic, one would obtain the number in the lower right-hand corner. In 
view of the many uncertainties in the analysis of the experimental data, it is impossible 
to tell whether the discrepancy between simulation and experiment is significant. An 
estimate of 7 based on the experimental enthalpy of fusion following Turnbull, would 
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yield 7 ~ llb[erg cm~'^] (see e.g. [H]). 

Huang et al. [9J estimated the solid-liquid surface free-energy density of NaCl from 
the nucleation rate at 550K. To achieve this, Huang et al. assumed that the CNT 
expression for the nucleation rate is valid. Under those assumptions, they obtained: 
7=119. 6[er(y' cm~^] for a spherical nucleus. 

In the present work, we can compute absolute nucleation rates without making 
use of CNT. The only assumption we make is that the barrier crossing is diffusive and 
that the Zeldovitch (pre)factor is well approximated by the form given in eqn.(jll). The 
Zeldovitch factors were found to be respectively Zi=0.016 for Ti and ^2=0. 013 for T2. 
The true Zeldovitch factor may be slightly different, but is in any event expected to be of 
(9(10~^). From our MD simulations we obtained the following estimates for the forward 
rates (eq. [71) : /^ji=0.013 ps~^ for Ti and /^jj=0.033 ps~^ for T2. Combining this 
information, we can compute the kinetic prefactor of eq. El k{Ti)= 6.9x10^^ m~^s~^ 
and k(T2)= 1.5x10^^ m~^s~^. As is to be expected, the kinetic prefactor depends only 
weakly on temperature. Using eqn.(|lj) we then calculated the nucleation rate. The 
results are: /(Ti)=3x lO^^^i m'^s'^ and /(T2)=4x 10^^^^ m-^s'^ 

These nucleation rates are about ten orders of magnitude higher than the estimated 
experimental rate at 905K {O{10^^)m~^s~^). Such a difference is hardly surprising 
because the nucleation rate is expected to increase rapidly with increasing supercooling. 

We also computed the nucleation rate using an algorithm based on the path- 
sampling techniques of refs. fldi [16] (Forward-Flux Sampling (FFS) and Transition- 
Interface Sampling (TIS)) . The value obtained at Ti = 800K is Iffs-tis(Ti) = 
0(10^"^^^) m~^s~^, which agrees surprisingly well with the one obtained using the 
diffusive barrier-crossing approach. 

As the path-sampling method does not depend on the choice of the reaction 
coordinate and does not require prior knowledge of the phase space density, we can 
conclude that the method based on the free-energy calculation gives us a good estimate 
for the nucleation rate at this temperature. 

We have also computed the nucleation rate at T=750 K, using the FFS-TIS method. 
The computed nucleation rate is IpFS-TisiT = 740i^')=O(10^^) m~^s~^. This is of the 
same order of magnitude as the nucleation rate obtained calculated by Koishi et al. [10] 
using "brute-force" MD. We cannot use the diffusive barrier crossing method at this 
temperature, as the barrier is too low to avoid spontaneous nucleation during long runs. 
However, if we assume that the kinetic pre-factor and the surface free-energy density 
do not vary much with temperature, we can use CNT to extrapolate the nucleation 
rate from 800K to 740K. We find: hxtrapiT = 740K) = O{10^°) m^^-^ which is 
considerably lower than the results of the direct calculations. This suggests that an 
extrapolation procedure based on CNT is not reliable. A summary of our numerical 
results for the nucleation barriers and rates are given in Table [3] 

A similar problem occurs if we try to extrapolate our numerical data at SOOi^' and 
825K to 905K, the temperature of the experiments of ref. [6]. If we can extrapolate 
our simulation results to the experimental temperature of 905 K (/3A/i=0.3) we obtain 
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T[K] 


/3A/i 


/3AGcrit 


f+ 

J ncrit 


I[m ^] 


lFFS-Tis[m ^] 


800 


0.54 


24 


0.013 


3 xl026±i 




825 


0.48 


29 


0.033 


4X1024±1 





Table 3. Summary of the simulation results for the calculation of the free-energy 
barrier and the nucleation rate for Tosi-Fumi NaCl. 



an estimated nucleation at 905 K that is 0(6 x 10^^)m^'^s^^. This is some nearly five 
orders of magnitude less than the experimentally observed rate. 

The discrepancy between simulation and experiment can be due to several reasons, 
a) There might be an appreciable (but unspecified) error in the experimental estimates 
(e.g. due to residual heterogeneous nucleation). b) The estimated error in the computed 
melting temperature of the Tosi-Fumi model is ±20K [20]. Such an uncertainty again 
easily translates into a variation of the nucleation rate by several orders of magnitude, c) 
In view of the extreme sensitivity of nucleation rates to the details of the intermolecular 
potential (see, e.g. [28]), the Tosi-Fumi potential may be inadequate to model nucleation 
in NaCl, even though it can reproduce the static properties of the solid and liquid NaCl 
[29] . d) Finally, it it is not quite correct to assume that the kinetic prefactor, the surface 
free energy and the latent heat of fusion are temperature- independent. 

We can also compare our calculated nucleation rate at 800K to the rate estimated 
with CNT. In order to do that, we need to compute the kinetic pre-factor k,cnt that, 
using the CNT approximations, is [5]: 

f^CNT — ^ Pliq r5 • \o) 



The attachment rate of particles to the critical nucleus (f^cnt) takes into account the 
number of available attachment sites on the surface of a spherical nucleus (n^/^^j) and 
depends on the jump frequency for bulk diffusion (-D5/ A^), where A is the atomic 
jump distance. Since the functional form of the nucleation barrier can be fitted to 
the corresponding CNT expression, the computed Zeldovitch factor (Z) coincides with 
the predicted one. We computed the self-diffusion coefficient with MD simulations 
using the DL_POLY package |23j in the supercooled liquid at Ti=800K. We found 
Dg"- = 3.4 X 10~^cm^s~^, in good agreement with an estimate {Dg°- = 2.3 x 10~^cm^s~^) 
based on extrapolation of the available experimental data of Ref. [30] to the temperature 
Ti. Since Dg°' /D^^ ^ 1, we only considered the self-diffusion of the Na"*" ions. We 
estimated A as a fitting parameter from the f^cru previously calculated; we obtained 
A(Ti) = 10"^ A. However, considering that the ion size is a^a ~ 1-1 A, this value for the 
jump distance seems unphysical (A ~ lOOo"). Typically, one would expect A to be of the 
order of a mean free path. In a molten salt, the mean free path of an ion is certainly 
less than a particle diameter. This discrepancy also suggests that the CNT picture is 
inadequate to describe crystal nucleation of NaCl. 

In summary, we have computed the crystal nucleation rate of sodium chloride 
from the melt using two independent methods: one based on calculations of the free- 
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energy barrier and the diffusive barrier-crossing and another based on a path-samphng 
approach. We have found that, to within an order of magnitude, the two approaches 
yield the same value for the nucleation rate. When we use Classical Nucleation 
Theory to extrapolate our numerical data to lower temperatures, we observe serious 
discrepancies with the results of direct calculations. When we use CNT to extrapolate 
to high temperatures, we find serious discrepancies with the nucleation rates found in 
experiments. Several factors may contribute to this discrepancy but, at present, it is 
not yet known which factor is most important. 
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Appendix A. Identification of crystalline clusters 

To distinguish between solid-like and liquid-like particles and identify the particles 
belonging to a solid cluster, we used the local bond order parameter introduced by 
Ten Wolde el al. [3l]. Although the method we used is the same as the one proposed 
by Ten Wolde, the definition of a 'solid-like' particle is not rigorously the same. The q 
vector and the thresolds selected were optimazed for the NaCl Tosi-Fumi model. 

First we computed a normalized complex vector ^4 for every particle i. Each 
component of this vector was given by: 

1 -^Nf^ii) 



Where Nb{i) is the number of neighbours of the particle i within a cut-off radious of AA 
(the first minimum in the Na-Cl radial distribution function). 

Then we computed a scalar product 54 (z) ■ ql{j) for every particle i with each of its 
neighbours particle j. A particle was considered to be 'solid-like' when at least 6 of the 
scalar products were bigger than 0.35. Finally two 'solid-like' particles were considered 
to be neighbours in the same cluster if they were closer than 3.4 A. 

Ionic fiuids are more ordered than Lennard- Jones or Hard Spheres ones, as the 
radial distribution function shows. However, with the method implemented by Ten 
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Wolde, we were able to clearly distinguish between solid-like and liquid-like particles in 
the NaCl Tosi-Fumi model. We enclose a plo HAll that shows the distributions of the 
number of scalar products bigger than 0.35 at T=800 K for the solid and for the liquid: 
for values bigger than 6 a particle was considered to be solid-like. 

Distributions of tlie number of conections at 800 K 




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

n 

Figure Al. Distributions of the number of scalar products bigger than 0.35 at T=800 
K for the sohd and for the hquid. 
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